"""This simulation shows for system dot(x)=1/(x-1)+g(t) with initial state x(0)>1
the system will never hit x=1, for even nonsmooth g(t)
"""

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np
 
def Alpha(d,mu,alpha_max):
    def f(s):
        if s>mu:
            o=0
        elif s<d:
            o=alpha_max
        else:
            o=1/(s-d)-1/(mu-d)
            if o>alpha_max:
                o=alpha_max
        return o
    return f


def f(t, y,alpha):
    x=y[0]
    v=y[1]

    phi=alpha(x)
    u=phi-10*(np.sign(t%10-5)+1)
    return [v, u]

if __name__ == '__main__':
    d=1;mu=3
    alpha=Alpha(d,mu,1e4)
    plt.figure()
    s=np.linspace(0,10,100)
    plt.plot(s,np.fromiter([ alpha(i) for i in s],dtype=float))
    plt.semilogy()
    plt.grid()
    plt.title(r'$\alpha(s)$')


    sol = solve_ivp(f, [0, 10], [2, 0], t_eval=np.linspace(0, 10, 10000),args=(alpha,))
    plt.figure()
    plt.subplot(211)
    plt.plot(sol.t, sol.y[0])
    plt.title(r'$x(t)$')
    plt.grid()
    plt.subplot(212)
    plt.plot(sol.t, 1/(sol.y[0]-1))
    plt.semilogy()
    plt.title(r'$1/(x(t)-1)$')
    plt.show()
